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SUMMARY 


A microcomputer system being developed by the authors is introduced. The 
parallel effort of compiling a series of compact finite element analysis pro- 
grams enables the execution of most computation on inexpensive microcomputers. 
The system is practically maintenance free and can be sustained by individual 
laboratories of standard scale in the educational or academic environment. 

As for the programs, *FEMN is discussed in some detail. The program is an 
extended version of the original linear analysis program FEM4 and is being 
tested for application to problems with material nonlinearities. 


INTRODUCTION 


The finite element analysis has reached the stage where the execution 
of the structural analysis is often considered routine. This is the case 
particularly in the industrial environment. However, so far the execution 
has relied largely on expensive hardware or costly remote time-sharing 
services. The role of the giant main-frame or super-computer in the solu- 
tion of large scale problems, e.g. the inelastic analysis of pressure vessels 
and piping systems operated at elevated temperature, will not be changed 
even in future applications. But it has been an ambition of engineers to 
perform a great portion of their analysis jobs on inexpensive and hopefully 
personal computers and thus be freed from being slave to the large systems. 
The development of microcomputer and associated finite element analysis pro- 
grams is a breakthrough in realizing this goal. 

The microcomputer system should be stand alone and almost maintenance 
free so that it can be sustained by individual laboratories especially in 
the educational or academic environment. The medium-sized engineering 
problems should be solved within a reasonable time limit and the system 
could also be adapted to multi-purpose usages, i.e. interactive compilation 
of fundamental computation routines, data management, preparation of engi- 
neering documents and reports, letter writing and so forth. In the present 
paper a compact system is introduced which is being built by the authors. 

In a parallel effort, a series of microcomputer finite element analysis 


277 



programs are being developed. The original version is FEM4 which is an 
elastic analysis program of plane stress, plane strain and axisymmetric 
problems (ref. 1, 2). It is extended to the nonlinear analysis program FEMN 
by an addition of the restart capability. The results of this innovation 
are manifold. By an incorporation of user defined subroutine MTKLN speci- 
fying material data, problems with material nonlinearities can be easily 
handled. Mesh division can be modified in the course of computation and 
thus the simulation and/or pursuit, for. example, of crack development in 
fracture mechanics becomes easier. In the following, some details of 
the microcomputer system and the program organization of FEMN are described 
with the example solution of a simple pilot problem. 


MICROCOMPUTER SYSTEM STRUCTURE 


Figure 1 illustrates the structure of the microcomputer system almost 
completed at the time. May 1980, of writing this paper. Zilog Z80 is used 
as the 8 bit CPU (Central Processing Unit) , and the capacity of main memory, 
which is composed of a ROM (Read Only Memory) and several RAM (Random Access 
Memory) boards, totals 64 kilobytes. The transfer of control, address and 
data between the components of the system is performed exclusively via S-lOO 
bus. For the purpose of connection and communication a number of interfaces 
are installed. The 8" floppy disk drive constitutes the secondary memory for 
mass storage and provides the housing of a compound of operating system, sup- 
porting language, finite element analysis and other computer programs. The 
standard disk operating system CP/M is used so that the problems in software 
exchange can be avoided. At the moment, program languages are BASIC and 
FORTRAN. It should be noted that two 240K dual disk systems are combined for 
commanding four floppy disk assemblage, although a dual system kit suffices 
to perform the standard operation. The authors intend to shorten the overall 
analysis time by an adoption of parallel processing that uses several CPUs 
and disk drives. The contemplated inclusion of the hard disk will increase 
the capacity of secondary memory to a great extent and may open the way to a 
novel system based on 16 bit microprocessors. 

Among the peripherals for I/O (Input /Output) purposes shown in figure 1, 
CRT unit and printer are the essentials- The function of CRT unit is manifold, 
as it is used for input of the data, interactive operation of the system, com- 
pilation of program segments and/or. subroutines, temporary display of the com- 
puted results, preparation of documents, e.g. the users’ manual, and so on. 

The prepared data, the completed list of the programs and the results of com- 
putation can be plugged into the printer for permanent recording. Plotter 
and graphics terminal are optional, but both are useful for the finite 
element modeling and post-processing graphics, e.g. automatic mesh genera- 
tion, model editing and plotting the computed results. 

The general purpose programs compiled to date for mounting on the micro- 
computer system are COMPOL, COMPOS, CALM, FEM4, MICRO-FEM and FEMN. The 
first two, written in BASIC, are essentially the microcomputer version of 
’COMPOsite material computation’ being developed by Tsai on a magnetic card 
calculator (ref. 3, 4). In the program names, L and S stand for the laminate 
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and sandwich composite structure respectively. CALM is a matrix operation 
program which is basically an interactive version of the first group opera-- 
tions in the program CAL (ref. 5, 6). FEM4 and its microcomputer version 
MICRO-FEM are prepared to solve the plane stress, plane strain and axisym- 
metric problems and then converted to FEMN to conduct nonlinear analysis. 
Restart capability is implemented so that the inelastic material behavior 
can be handled. Complementary modification is an addition of the user 
defined subroutine MTRLN specifying the material properties which were for- 
merly input through element data cards. The users who tailor the subroutine 
MTRLN according to their material data can perform conveniently the inelastic 
analysis on the microcomputer. An example of MTRLN is shown in the following 
section. 


ORGANIZATION OF NONLINEAR ANALYSIS PROGRAM FEMN 


FEMN is composed of two parts FEMAB and FEMCD which are concerned with 
the preparation of input data file and the solution procedure of the prob- 
lem. The major feature of the program is that it uses d 3 mamic storage 
allocation which means the complete elimination of common statement. This 
function is performed by subroutines OPENS, CLOSE, PSEEK and POOLWT as shown 
in figures 2a and 2b. 

The program organization of FEMAB is shown in figure 2a. It reads title 
and control cards first. Then follows the input of node and element data from 
which the index or integer joint array is formed and stored on IFIL2, IFIL2 
accommodates also load data and IFIL6 is essentially a storage of input and 
processed element data. The formation of strain-displacement matrix iB from 
the data in IFIL6 is performed by subroutine MGN and the result is written 
on IFIL3 to be read in FEMCD subsequently. Finally the initial displacement 
data (usually the cleared zero displacement) are stored on IFIL5 for subse- 
quent updating by the solution obtained through FEMCD. In the following the 
principal functions of individual subroutines in FEMAB are summarized. 

PINPG Preparatory segment for subroutine INPUTG 

INPUTG Input generation, read input data in sequence and compile index or in- 
teger joint array 
RNODE Read node data 

RELEM Read element data including material specification number and element 
thickness 

MKIDX Make index for attributing merging point and coordinate to input and 

processed element data and also create index for skyline assemblage of 
stiffness matrix 

RLOAD Read loading step sequence and nodal force and/or displacement data 

OPENS Open storage area for array from bottom of POOL in main memory 
CLOSE Close a part or whole storage area by deleting unused array 
PSEEK Search for array data by its name 

POOLWT Debug write wanted array data in POOL for inspection 

PMGN Preparatory segment for subroutine MGN 

MGN Command sequential generation of element matrices 

ISOBMN Generate element strain-nodal displacement matrix B for 3-8 variable 
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STDMA 


node parametric element 

Evaluate components of strain-nodal displacement matrix B and determinant 
of associated Jacobi matrix J 
DISPI Initialize displacements, by either clearing for initial run or 

entering displacement data of preceding computation for restart run 

After the execution of FEMAB described above, the next program block FEMCD 
is called by the main of FEMN, Fxmctions of FEMCD, whose organization is shown 
in figure 2b, are the formation and assemblage of the element stiffness matrix 
and subsequently the solution of the problem. Material data are input via the 
user defined subroutine MTRLN and then the stress-strain matrix D for each 
incremental stage of loading is evaluated in subroutine DMXMKN. FEMCD starts 
its operation by a transfer or reading of the data stored on IFIL5, which are 
title of the problem, integer data, initial cleared displacement and load data 
or their values obtained in preceding step of loading sequence. The stress- 
nodal displacement matrix S and element stiffness matrix K are formed in sub- 
routine SMXMKN, the latter being stored in the appropriate locations in the 
overall stiffness matrix by referring to the index prepared on IFIL2. IFIL3 
and 4 are used as the seesaw external memory for integer data and the element 
strain-displacement matrix B. Skyline or profile active column method of data 
acquisition is used for saving area in the main memory. Therefore all sub- 
routines prefixed by capitals SK in SOLVEN, SKDCNP etc., take advantage of 
the sky lined form of storage for the manipulation of data. Newton-Raphson 
iteration procedure is incorporated in subroutine SOLVEN, some details of which 
are discussed in the next section concerned with the solution of an elementary 
sample problem. The following summarises the function of individual sub- 
routines relevant to FEMCD. 

PSOLVN Preparatory segment for subroutine SOLVEN 

SOLVEN Solve overall stiffness equation for unknown nodal displacement and com- 
pute reaction at constrained node; iterative procedure is incorporated 
in this solver for nonlinear problems 
VECTWN Print out computed displacement and reaction vector at each stage of 
loading 

SMXMKN Evaluate stress-nodal displacement matrix S and element stiffness 
matrix K, synthesize overall stiffness matrix by referring to 
merging point index stored on IFIL2, and also determine equivalent 
nodal force from current stress data for equilibrium check 
SWRITE Write components for debugging of active columns in matrix S stored 
by skyline method 

BWRITE Write components of vector B for debugging 

SKDCNP Cholesky decomposition of symmetric positive definite matrix by skyline 
method 

SKXMLU Multiply, add and/or subtract matrix components in skyline storage 
CONVCK Check convergence of solution being obtained by Newton-Raphson iterative 
procedure 

SKFWD Forward elimination by skyline method 
SKBKW Backward substitution by skyline method 

DMXMKN Evaluate components of stress-strain matrix D of constitutive equation 
MTRLN User defined subroutine specifying elastic and inelastic material prop- 
erties 
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STRSUM Add stress/strain increments to update values of stress/straln 
PRINST Compute principal stresses and their directions 

STRPRN Print out stress/strain solutions at respective Gauss integration 
points, together with coordinates of Gauss points 

Finite element used in FEM4 as well as in FEMN is 4-8 variable nodes para- 
metric quadrilateral with the following interpolation functions (ref. 7). 

For corner nodes 1-4 

Ni = Si - (Ng+N^)/2, - (N^+Ng)/2 

Ng - Sg - (Ng+N^)/2, - (N^+Ng)/2 

for midedge nodes 5-8 



where stands for the trial functions defined as (ref, 8, 9) 

= (l-C)(l-n)/4, = (l+?)(l-n)/4 

Sg = (l+?)(l+n)/4, = (l-?)(l+n)/4 

= (l-52)(l-n)/2, Sg = (l+5)(l-n2)/2 

= (l-? 2 )(i+n)/ 2 , Sg = (l-5)(l-n2)/2 

By coalescing an edge of the quadrilateral to a single point a triangular ele- 
ment is produced. It can be shown (ref .8, 9) that the resulting element coin- 
cides with the conventional constant stress/strain element when the primary 
quadrilateral element is four-noded. The number of integration points in 
Gauss quadrature can be one to five by five in accordance with the users’ 
specification. 

The input card or data sequence in FEMN is summarized in table I. An 
example of input data preparation as well as the user defined subroutine MTRLN 
is illustrated in the next section. 

SOLUTION OF SAMPLE NONLINEAR PROBLEM AND REMARKS 


As a sample problem, nonlinear behavior of a composite block specimen 
shown in figure 3 is analysed. The block consists of an ideally plastic ele- 
ment 101 and a nonlinear one 501 with a negative slope segment at a large 
strain as depicted in figure 4 based on the material data of figure 3. Loading 
sequences are summarized in figure 3 and the solid curve in figure 4 is the 
theoretical load-displacement relation of the block under axial tensile load- 
ing. It is noted that the loading condition in numerical analysis is given by 
the force increment for step 1-3 and 7-9, while in step 4-6 it is given by the 
displacement increment . 

Table II is the image of input cards prepared for the solution of the 
sample problem and serves to illustrate simplicity of the data preparation. 

281 



Specifically table II is concerned with the first loading sequence, l.e. 
step 1-3. The solution for consecutive loading conditions, step 4-6 and 7-9 
in the present example, is obtained by restarting the execution with the 
renewal of input data and the use of the solution obtained in the preceding 
step and stored on an appropriate file. Table III is the subroutine MTRLN 
written for this sample problem. The program FEMN is versatile because the 
user can easily tailor the subroutine MTENL so that it characterizes particu- 
lar nonlinear properties of the material of interest. It must be emphasized 
that the anisotropic material behaviors are easily incorporated in the program. 

Figure 5 depicts solution convergence in the sample problem. The itera- 
tive procedure that the present version of FEMN employs is a modified Newton- 
Raphson method with incorporation of equivalent nodal force {F^}. It com- 
pensates the imbalance of force equilibrium at the nodes and is given by 

{F^} = {F} - |[B]'^{a}dV 

where {F} denotes the prescribed nodal force, [B] and {a} are the strain-nodal 
displacement matrix and the current stress. Convergence is satisfactory in the 
present example and it should be noticed that the computed results lie on the 
theoretical curve exhibiting sharp turning points. 

Test of convergence in case of the large scale problem, sophistication of 
iterative procedure and extension of the program to three dimensions are the 
next steps that are to be taken. Moreover, the development of parallel pro- 
cessing and the installation of suitable hard disk will increase the speed and 
capacity of the system. 
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TABLE I 


INPUT DATA FORMAT OF FEMN 


(1) TITLE CARD (18A4) 

COL 1-72 PROBLEM IDENTIFICATION ETC BY ALPHANUMERIC CHARACTER 

(2) CONTROL CARD (2311) 

COL 1 =0 AXISYMMETRIC 

=1 PLANE STRAIN 

=2 PLANE STRESS 

2 NUMBER OF INTEGRATION POINTS (1-5) FOR GAUSSIAN QUADRATURE 

3* =0 INITIAL START 

>0 RESTART 

4- 5* NUMBER OF ITERATION IN NEWTION-RAPHSON METHOD 

6- 20 BLANK 

21 >0 DEBUG WRITE IN MODULE INPUG 

22 >0 DEBUG WRITE IN MODULE MG 

23 >0 DEBUG WRITE IN MODULE SOLVEN 

(3) NODE HEADER CARD (A4) 

COL 1-4 'NODE' 

(4) NODE DATA CARDS (A4, 16, 2F10.0, lOX, FIO.O, lOX, 211) 

COL 1- 4 'NODE' 

7- 10 NODE NUMBER 
11-20 X(R) COORDINATE 
21-30 Y(Z) COORDINATE 

41-50 OBLIQUE ANGLE (DEG) OF LOCAL COORDINATE 

61 =1 X-DOF CONSTRAINED OR X(R)-DISPL GIVEN 
=0 OR BLANK FREE OR X(R)-LOAD GIVEN 

62 =1 Y-DOF CONSTRAINED OR Y(Z)-DISPL GIVEN 

=0 OR BLANK FREE OR Y(Z)-LOAD GIVEN 

(5) ELEM HEADER CARD (A4) 

COL 1- 4 'FLEM' 

(6) ELEM DATA CARDS (A4, 16, 815, 15*, 215*, F8.0) 

COL 1- 4 'ELEM' 

7-10 ELEMENT NUMBER 
11-15 1ST NODE NO. 

16-20 2ND NODE NO. 

21-25 3RD NODE NO. 

26-30 4TH NODE NO. 

31-35 5TH NODE NO. 

36-40 6TH NODE NO. 

41-45 7TH NODE NO. 

46-50 8TH NODE NO. 

51-55* MATERIAL SPECIFICATION NUMBER 
56-65* FOR EXTENSION OF PROGRAM BY USERS 
66-73 ELEMENT THICKNESS 


- 
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(7) LOAD HEADER CARD (A4) 

COL 1- 4 'LOAD' 

(8) LOAD OR DISPLACEMENT STEP CARD (A4, 16, 6F10.0)* 

COL 1- 4 'STEP' 

5-10 STEP NUMBER 

11-70 FOR PROGRAM EXTENSION 

(9) LOAD OR DISPLACEMENT DATA CARD (A4, 16, 2F10.0) 

COL 1- 4 'LOAD' 

5-10 NODE NUMBER ON WHICH GIVEN LOAD OR DISPLACEMENT IS APPLIED 
11-20 X(R) GIVEN NODAL FORCE OR DISPLACEMENT 

21-30 Y(Z) GIVEN NODAL FORCE OR DISPLACEMENT 

(10) END CARD (A4) 

COL 1- 3 'END' 

* INDICATES ADDITION OR MODIFICATION APPLIED TO LINEAR ANALYSIS PROGRAM FEM4 
AND/OR MICRO-FEM 


TABLE II INPUT DATA IMAGE OF SAMPLE PROBLEM OF FIGURE 3 


2-ELEMENT NONLINEAR MODEL TEST 1980-5-16 (ITER MAX= 8) 


22 8 




NODE 




NODE 

1 

0.000 

0. 000 

NODE 

2 

0. 000 

20.000 

NODE 

3 

0.000 

40. 000 

NODE 

11 

25.000 

0.000 

NODE 

13 

25. 000 

40.000 

NODE 

21 

50. 000 

0. 000 

NODE 

22 

50.000 

20.000 

NODE 

23 

50.000 

40.000 

ELEM 




ELEM 

101 

1 21 

23 3 

ELEM 

501 

1 21 

23 3 

LOAD 




STEP 

1 



LOAD 

21 

4000.000 

0. 000 

LOAD 

22 

16000.000 

0.000 

LOAD 

23 

4000. 000 

0. 000 

STEP 

2 



LOAD 

21 

2000. 000 

0.000 

LOAD 

22 

8000. 000 

0.000 

LOAD 

23 

2000. 000 

0. 000 

STEP 

3 



LOAD 

21 

2000.000 

0. 000 

LOAD 

22 

8000.000 

0. 000 

LOAD 

23 

2000. 000 

0.000 

END 





11 

10 

10 

01 

00 

01 

00 

00 

11 22 13 2 0 0. . 

11 22 13 2 1 0 . 


(65) 

— > 0 6.000 
0 6.000 
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n n n n n n a 


TABLE III EXAMPLE OF USER DEFINED SUBROUTINE MTRLN 


SUBROUTINE MTRLN(KK, IK, EK, ITER, ISTP) 


. IK(4) MflTERIfiL IDENTIFICATION NUMBER 

. EK(IK(3:0 COMPUTED STRAIN VALUE 
. EKaK(B)) CORRESPONDING STRESS VALUE 
. ITER=0, ROUTINE DETERMINES TANGENT MODULUS 
. ITER>0, ITERATES STRESS FOR COMPUTED STRAIN VALUE 


DIMENSION KKd), IK<l),EKa) 

MID=IK<4:) 

IKS=IK(S> 

IKS=IK(3:) 

IFUTER. GT. 0) GO TO GOOO 
C +.+.-*•+ SLOPE OF STRESS-STRAIN CURVE 

IF(MID.GT.O) GO TO 3000 
C PERFECTLY PLASTIC 

EKC2)=0. 3 

IF(EK(IK3). GT. 10. OE-3:) GO TO 2100 
EK(1 )=2. 0E4 
GO TO 3000 
2100 EK<1>=0.0 
GO TO 3000 

C NONLINEAR MATERIAL + + 

3000 EK<2>=0.3 

IFCEKdKO) . GT. 3. 0E-3> GO TO 3100 
EK<1)=2. 0E4 
GO TO 3000 

3100 IFCEKCIKO) . GT. 6. 0E-3) GO TO 3200 
EK< 1 >=0. 0 
GO TO 3000 
3200 EK(1>=-1. OE+4 
GO TO 3000 

C + STRESS VALUE FOR COMPUTED STRAIN 4=*-** 

GOOO IFCMID. GT. 0) GO TO 7000 
C .+.+-+:4: PERFECTLY PLASTIC 

IFCEKdKS:) . GT. 10. OE-3) GO TO G200 
EKdKS)=2, 0E4+EKdK3) 

GO TO 3000 
G200 EK (IKS) =200.0 
GO TO 3000 

C *+-.+:+: NONLINEAR MATERIAL 

7000 IF(EKdK3) . GT. 3. OE-3) GO TO 7200 
EKdKS)=2. 0E4+EKdK3) 

GO TO 3000 

7200 IF(EKdK3) . GT. G. OE-3) GO TO 7400 
EKdKS)=GO. 0 
GO TO 3000 

7400 EK (IKS) =120. 0-1 . 0E4*EK (IK3) 

3000 RETURN 
END 
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Figure 1.- Structure of microcomputer system. 
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(a) FEMAB. 



(b) FEMCD. 

Organization of FEMAB and FEMCD. 
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